Exercises Laboratory Session 05

knitr::opts_chunk$set(dev = "ragg_png")
options(digits=5) # set number of digits equal to 5

Packages and functions

# tidyverse
library(tidyverse)

# others
library(gridExtra)
library(R6)
library(kableExtra)
library(glue)
library(patchwork)
library(latex2exp)
library(metR)   # label contour plot

Exercise 1

The number of particles emitted by a radioactive source during a fixed interval of time (\(\Delta t=10s\)) follows a Poisson distribution on the parameter \(\mu\). The number of particles observed during consecutive time intervals is: 4, 1, 3, 1 and 3.

a) Suppose a uniform prior distribution for the parameter \(\mu\)

Determine and draw the posterior distribution for \(\mu\), given the data

Assuming a uniform prior, the posterior is simply proportional to the likelihood. Being a Poisson process, the posterior then is a \(Gamma(\alpha, \lambda)\), where \(\alpha = \sum_{j=1}^{n} y_j+1\) and \(\lambda=n\).

mu <- seq(0, 10, length=1000)

n.particles <- c(4, 1, 3, 1, 3)
n <- length(n.particles)

unif.shape <- sum(n.particles+1)
unif.rate <- n
unif.posterior <- dgamma(mu, shape=unif.shape , rate=unif.rate)

ggplot()+
  geom_line(aes(x=mu, y=unif.posterior), color='steelblue', size=0.7)+
  labs(
    title = 'Posterior distribution, positive uniform prior',
    x = TeX('$\\mu$'), 
    y = TeX('$P(\\mu|\\{y_i\\})$')
  )+
  theme_bw()

Evaluate mean, median and variance, both analytically and numerically in R

The median can not be computed analytically, as there is not a simple closed form expression for it.

argmax <- function(x){
    max.x <- max(x)
    found <- FALSE
    
    i <- 1
    while (!found & i<=length(x)){
      if (x[i]==max.x){
        found <- TRUE
        return(i)
      }
      i <- i+1
    }
  }

gamma.mean.analytically <- function(shape, rate){
  return(shape/rate)
}
gamma.variance.analytically <- function(shape, rate){
  return(shape/rate^2)
}
gamma.mode.analytically <- function(shape, rate){
  return((shape-1)/rate)
}
gamma.mode.numerically <- function(param, distr.data){
  return(param[argmax(distr.data)])
}
gamma.median.numerically <- function(shape, rate){
  return(qgamma(0.5, shape=shape, rate=rate))
}
gamma.mean.numerically <- function(shape, rate, upp.bound=30){
  return(integrate(function(x){dgamma(x, shape=shape, rate=rate)*x}, 0, upp.bound)$value)
}
gamma.variance.numerically <- function(shape, rate, upp.bound=30){
  k2 <- integrate(function(x){dgamma(x, shape=shape, rate=rate)*(x^2)}, 0, upp.bound)$value
  k1 <- gamma.mean.numerically(shape, rate, upp.bound)
  
  return(k2-k1^2)
}

glue('Median, computed numerically: mu={format(gamma.median.numerically(unif.shape, unif.rate), digits=4)}')

## Median, computed numerically: mu=3.334

df <- data.frame(row.names=c('Analytically', 'Numerically'))

df['Mean'] = c(gamma.mean.analytically(unif.shape, unif.rate), gamma.mean.numerically(unif.shape, unif.rate))
df['Variance'] = c(gamma.variance.analytically(unif.shape, unif.rate), gamma.variance.numerically(unif.shape, unif.rate))
df['Mode'] = c(gamma.mode.analytically(unif.shape, unif.rate), gamma.mode.numerically(mu, unif.posterior))

df %>% 
  kable(digits=3, caption='Posterior statistics, with uniform prior') %>%
  kable_classic(full_width=FALSE)

Posterior statistics, with uniform prior

	Mean	Variance	Mode
Analytically	3.4	0.68	3.200
Numerically	3.4	0.68	3.203

unif.mean <- gamma.mean.analytically(unif.shape, unif.rate)
unif.mode <- gamma.mode.analytically(unif.shape, unif.rate)
unif.std <- sqrt(gamma.variance.analytically(unif.shape, unif.rate))

b) Suppose a Jeffrey’s prior for the parameter \(\mu\)

Determine and draw the posterior distribution for \(\mu\), given the data

jeffrey.shape <- sum(n.particles+1/2)
jeffrey.rate <- n
jeffrey.posterior <- dgamma(mu, shape=jeffrey.shape, rate=jeffrey.rate)

ggplot()+
  geom_line(aes(x=mu, y=jeffrey.posterior), color='steelblue', size=0.7)+
  labs(
    title = 'Posterior distribution with Jeffrey’s prior',
    x = TeX('$\\mu$'), 
    y = TeX('$P(\\mu|\\{y_i\\})$')
  )+
  theme_bw()

Evaluate mean, median and variance, both analytically and numerically in R

glue('Median, computed numerically: mu={format(gamma.median.numerically(unif.shape, unif.rate), digits=4)}')

## Median, computed numerically: mu=3.334

df <- data.frame(row.names=c('Analytically', 'Numerically'))

df['Mean'] = c(gamma.mean.analytically(jeffrey.shape, jeffrey.rate), gamma.mean.numerically(jeffrey.shape, jeffrey.rate))
df['Variance'] = c(gamma.variance.analytically(jeffrey.shape, jeffrey.rate), gamma.variance.numerically(jeffrey.shape, jeffrey.rate))
df['Mode'] = c(gamma.mode.analytically(jeffrey.shape, jeffrey.rate), gamma.mode.numerically(mu, jeffrey.posterior))

df %>% 
  kable(digits=3, caption="Posterior statistics, with Jeffrey's prior") %>%
  kable_classic(full_width=FALSE)

Posterior statistics, with Jeffrey’s prior

	Mean	Variance	Mode
Analytically	2.9	0.58	2.700
Numerically	2.9	0.58	2.703

jeffrey.mean <- gamma.mean.analytically(jeffrey.shape, jeffrey.rate)
jeffrey.mode <- gamma.mode.analytically(jeffrey.shape, jeffrey.rate)
jeffrey.std <- sqrt(gamma.variance.analytically(jeffrey.shape, jeffrey.rate))

c) Evaluate a 95% credibility interval for the results obtained with both priors

Compare the result with that obtained using a normal approximation for the posterior distribution, with the same mean and standard deviation.

unif.cred.int <- c(qgamma(0.025, unif.shape, unif.rate), qgamma(0.975, unif.shape, unif.rate))
jeffrey.cred.int <- c(qgamma(0.025, jeffrey.shape, jeffrey.rate), qgamma(0.975, jeffrey.shape, jeffrey.rate))

gg_post_complete <- function(post.funct, limits, mode, Title='Posterior'){
  mu <- seq(0, 10, length=1000)
  mu95 <- seq(limits[1], limits[2], length=1000)
    
  ggplot()+
    geom_line(aes(x=mu, y=post.funct(mu)), size=0.7, color='steelblue')+
    geom_area(aes(x=mu95, y=post.funct(mu95)), fill='lightblue', color='steelblue', size=0.7)+
    labs(
      title=Title, 
      x=TeX('$\\mu$'), 
      y=TeX('$P(\\mu|\\{y_i\\})$')
    )+
    ylim(0, 0.6)+
    geom_vline(xintercept=limits, linetype='dashed', size=0.6)+
    geom_vline(xintercept=mode, linetype='dashed', size=0.6, color='firebrick')+
    annotate('text', x=limits[1]-0.6, y=0.55, label=paste('x1=', format(limits[1], digits=2, nsmall=2), sep=''))+
    annotate('text', x=limits[2]+0.6, y=0.55, label=paste('x2=', format(limits[2], digits=2, nsmall=2), sep=''))+
    annotate('text', x=mode+0.4, y = 0.58, label=paste(format(round(mode, 2), nsmall = 2), sep=''), color='firebrick')+
    theme_bw()
}

df <- data.frame(row.names=c('Uniform prior', "Jeffrey's prior"))
df['Lower bound']=c(unif.cred.int[1], jeffrey.cred.int[1])
df['Upper bound']=c(unif.cred.int[2], jeffrey.cred.int[2])

df %>%
  kable(caption='95% credibility interval', digits=3) %>%
  kable_classic("striped", full_width = F)

95% credibility interval

	Lower bound	Upper bound
Uniform prior	1.981	5.197
Jeffrey’s prior	1.605	4.572

gg_post_complete(function(x){dgamma(x, shape=unif.shape , rate=unif.rate)}, unif.cred.int, unif.mode, Title='Posterior distribution, positive uniform prior')

gg_post_complete(function(x){dgamma(x, shape=jeffrey.shape , rate=jeffrey.rate)}, jeffrey.cred.int, jeffrey.mode, Title="Posterior distribution, Jeffrey's prior")

Using the normal approximation:

gg_post_complete(function(x){dnorm(x, mean=unif.mean, sd=unif.std)}, 
                 c(qnorm(0.025, mean=unif.mean, sd=unif.std), qnorm(0.975, mean=unif.mean, sd=unif.std)),
                 unif.mean, 
                 Title='Posterior distribution using the normal approx')+
  labs(subtitle='Mean and std from the posterior with the uniform prior')

gg_post_complete(function(x){dnorm(x, mean=jeffrey.mean, sd=jeffrey.std)}, 
                 c(qnorm(0.025, mean=jeffrey.mean, sd=jeffrey.std), qnorm(0.975, mean=jeffrey.mean, sd=jeffrey.std)),
                 jeffrey.mean, 
                 Title='Posterior distribution using the normal approx')+
  labs(subtitle="Mean and std from the posterior with the Jeffrey's prior")

Using the normal approximation we obtain very similar results. For example, with the positive uniform prior, the mode and the mean of the posterior distribution are, respectively, 3.2 and 3.4, while using the normal approximation, we obtain a mode and a mean of 3.4. Similarly, also the credibility intervals are almost the same. This is true also when considering the Jeffrey’s prior.

Notice that this results derive from the Central Limit Theorem, that in this case can be applied even if the number of measurements is quite small (5).

Exercise 2

Given the problem of the lighthouse discussed last week, study the case in which both the position along the shore (\(\alpha\)) and the distance out at sea (\(\beta\)) are unknown.

Exercise 3

Given the Signal over Background example discussed last week, analyze and discuss the following cases:

• a) Vary the sampling resolution used to generate the data, keeping the same sampling range xdat <- seq(from=-7*w, to=7*w, by=0.5*w)
  • change the resolution $w = {0.1, 0.25, 1, 2, 3}
  • check the effect on the results
• b) Change the ratio A/B used to simulate the data (keeping both positive in accordance with the prior)
  • check the effect on the results

Recall the signal shape is \[ S_k = \Delta t\left\{A\cdot\text{exp}\left[-\frac{(x_k-x_0)^2}{2w^2}\right]+B\right\} \] where \(\Delta t\) is the exposure time, \(B\) and \(A\) are the background and signal amplitudes, \(x_0\) is the centre of the signal peak and \(w\) its width.

Define general classes and methods

signal <- R6Class("signal",
    public = list(
      x  = NULL,
      A  = NULL,
      B  = NULL,
      x0 = NULL,
      w  = NULL,
      delta.t = NULL,
      signal = NULL, 
      ddat = NULL, 

      initialize = function(x=NA, A=NA, B=NA, x0=NA, w=NA, delta.t=NA) {
        self$x  <- x
        self$A  <- A
        self$B  <- B
        self$x0 <- x0
        self$w  <- w
        self$delta.t <- delta.t
        self$signal  <- self$delta.t * (self$A*exp(-(self$x-self$x0)^2/(2*self$w^2)) + self$B)
      }, 
      
      generate = function(){
        set.seed(205)
        self$ddat <- rpois(length(self$signal), self$signal)
      }, 
      
      plot.signal = function(){
        geom_line(aes(x=self$x, y=self$signal), col='navyblue', size=0.7)
      }, 
      
      plot.data = function(){
        if(is.null(self$ddat)) self$generate()
        
        if(self$A/self$B>=1.5) xdat.off <- self$x-(self$x[argmax(self$ddat)]+self$x[argmax(self$ddat)+1])/2
        else xdat.off <- self$x
        geom_step(aes(x=xdat.off, y=self$ddat), col='firebrick', size=0.7)
      }
    )
)

posterior <- R6Class("posterior",
    inherit = signal, 
    
    ## PUBLIC
    public = list(
      z  = NULL, 
      p_a_D  = NULL,
      p_b_D  = NULL, 
      p_a_bD = NULL, 
      p_b_aD = NULL, 
      
      log.unnorm.post.grid = function(){
        z <- matrix(data=NA, nrow=length(private$a.vect), ncol=length(private$b.vect))
        for(j in 1:length(private$a.vect)) {
          for(k in 1:length(private$b.vect)) {
            z[j,k] <- private$log.post(private$a.vect[j], private$b.vect[k])
          }
        }
        self$z <- z - max(z)
      },
      
      a.b.vectors = function(a.vect, b.vect){
        private$a.vect = a.vect
        private$b.vect = b.vect
      },
      
      plot.2d.posterior = function(title=''){
        if(is.null(self$ddat)) super$generate()
        
        self$log.unnorm.post.grid()
        
        df <- expand.grid(a=private$a.vect, b=private$b.vect)
        df$z <- as.numeric(exp(self$z))
        
        ggplot(df, aes(a, b, z=z))+
          geom_contour(bins=5, color='navyblue', size=0.7)+
          geom_text_contour(stroke = 0.2, rotate=F)+
          geom_vline(xintercept=self$A, col='firebrick', alpha=1, size=0.7)+
          geom_hline(yintercept=self$B, col='firebrick', alpha=1, size=0.7)+
          geom_vline(xintercept=0, col='black', linetype='dashed', size=0.7)+
          labs(
            x='amplitude, A', 
            y='background, B'
          )+
          {if (title!='') labs(title=title)}+
          xlim(c(self$A-1.5, self$A+0.5))+
          ylim(c(0.5, 1.5))+
          theme_bw()+
          {if (title!='') theme(plot.title  = element_text(size=18))}
      }, 
      
      marginalization = function(delta_a, delta_b){
        # Compute normalized marginalized posteriors , P(a|D) and P(b|D)
        self$p_a_D <- apply(exp(self$z), 1, sum) 
        self$p_a_D <- self$p_a_D/(delta_a*sum(self$p_a_D))
        self$p_b_D <- apply(exp(self$z), 2, sum) 
        self$p_b_D <- self$p_b_D/(delta_b*sum(self$p_b_D))
         
        # Compute normalized conditional posteriors , P(a|b,D) and P(b|a,D)
        self$p_a_bD <- exp(Vectorize(private$log.post , "a.i")(private$a.vect, self$B))
        self$p_a_bD <- self$p_a_bD/(delta_a*sum(self$p_a_bD))
        self$p_b_aD <- exp(Vectorize(private$log.post , "b.i")(self$A , private$b.vect))
        self$p_b_aD <- self$p_b_aD/(delta_b*sum(self$p_b_aD))
      }, 
      
      # Plot the 1D marginalized posteriors
      plot.marginalized = function(title=''){
        
        gg_marg_b <- ggplot()+
          geom_line(aes(x=private$b.vect, y=self$p_b_D), size=1, color='navyblue')+
          geom_line(aes(x=private$b.vect, y=self$p_b_aD), size=1, linetype='dashed')+
          geom_vline(xintercept=self$b, col='firebrick', size=1)+
          labs(
            x='background, B',
            y='P(B|D) and P(B|A,D)',
            title=title
          )+
          {if (title!='') labs(title=title)}+
          ylim(1.05*c(0, max(self$p_b_D, self$p_b_aD)))+
          theme_bw()+
          theme(legend.title = element_blank(), 
                axis.text   = element_text(size=16),
                axis.title  = element_text(size=16),
                plot.title  = element_text(size=18))

        gg_marg_a <- ggplot()+
          geom_line(aes(x=private$a.vect, y=self$p_a_D), size=1, color='navyblue')+
          geom_line(aes(x=private$a.vect, y=self$p_a_bD), size=1, linetype='dashed')+
          geom_vline(xintercept=self$a, col='firebrick', size=1)+
          labs(
            x='amplitude, A',
            y='P(A|D) and P(A|B,D)', 
            title=title
          )+
          {if (title!='') labs(title=title)}+
          ylim(1.05*c(0,max(self$p_a_D, self$p_a_bD)))+
          theme_bw()+
          theme(legend.title = element_blank(),
                axis.text   = element_text(size=16),
                axis.title  = element_text(size=16), 
                plot.title  = element_text(size=18))

        gg_marg_b+gg_marg_a
      }, 
      
      param_values = function(delta_a, delta_b, print_param=TRUE, return_param=FALSE){
        if(is.null(self$p_a_D)) self$marginalization(delta_a, delta_b)

        mean_a <- delta_a * sum(private$a.vect * self$p_a_D)
        mean_b <- delta_b * sum(private$b.vect * self$p_b_D)
        sd_a <- sqrt( delta_a * sum((private$a.vect-mean_a)^2 * self$p_a_D) )
        sd_b <- sqrt( delta_b * sum((private$b.vect-mean_b)^2 * self$p_b_D) )
        
        # Covariance normalization is performed with ’brute force ’. 
        # The normalization constant is Z = delta_a*delta_b*sum(exp(z)).
        # This is independent of (a,b) so can be calculated outside of the loops.
        cov_ab <- 0
        for(j in 1:length(private$a.vect)) {
          for(k in 1:length(private$b.vect)) {
            cov_ab <- cov_ab + (private$a.vect[j]-mean_a)*(private$b.vect[k]-mean_b)*exp(self$z[j,k])
        }
        }
        cov_ab <- cov_ab / sum(exp(self$z))
        rho_ab <- cov_ab / (sd_a * sd_b)
        
        if (print_param){
          cat("  a = ", mean_a, "+/-", sd_a, "\n")
          cat("  b = ", mean_b, "+/-", sd_b, "\n")
          cat("rho = ", rho_ab, "\n")
        }
        
        if (return_param){
          return(c(mean_a, mean_b, sd_a, sd_b, rho_ab))
        }
      }
    ), 
    
    ## PRIVATE
    private = list(
      a.vect = NULL, 
      b.vect = NULL, 
      
      log.post = function(a.i, b.i) {
        if(a.i<0 || b.i<0) {return(-Inf)} # effect of the prior
        else{
          sig.tmp <- signal$new(self$x, a.i, b.i, self$x0, self$w, self$delta.t)
          return(sum(dpois(self$ddat, lambda=sig.tmp$signal, log=TRUE)))
        
        }
      }
    )
)


gg_signal = function(sig1, sig2, title=''){
  # sig1 is the true signal to plot
  # sig2 contains the generated data
  ggplot()+
    sig1$plot.signal()+
    sig2$plot.data()+
    xlim(range(sig1$x))+
    ylim(range(c(sig1$signal, sig2$ddat)))+
    labs(
      x='x', 
      y='Signal+Background counts'
    )+
    {if (title!='') labs(title=title)}+
    theme_bw()+
    {if (title!='') theme(plot.title  = element_text(size=18))}
}

Reproduce the example done in class, to check if the above classes work correctly

x0      <- 0  # Signal peak
w       <- 1  # Signal width
A.true  <- 2  # Signal amplitude
B.true  <- 1  # Background amplitude
Delta.t <- 5  # Exposure time

sig.test <- posterior$new(seq(from=-7*w, to=7*w, by=0.5 *w), A.true, B.true, x0, w, Delta.t)
sig.plot <- posterior$new(seq(from=-7*w, to=7*w, by=0.05*w), A.true, B.true, x0, w, Delta.t)

gg_signal(sig.plot, sig.test)

# sampling grid for computing posterior
alim <- c(0.0, 4.0)
blim <- c(0.5, 1.5)
Nsamp <- 100
uniGrid <- seq(from=1/(2*Nsamp), to=1-1/(2*Nsamp), by=1/Nsamp)

delta_a <- diff(alim)/Nsamp
delta_b <- diff(blim)/Nsamp
a <- alim[1] + diff(alim)*uniGrid
b <- blim[1] + diff(blim)*uniGrid

# 2d posterior
sig.test$a.b.vectors(a, b)
sig.test$plot.2d.posterior()

sig.test$marginalization(delta_a, delta_b)
sig.test$plot.marginalized()

sig.test$param_values(delta_a, delta_b, print_param=TRUE)

##   a =  1.6302 +/- 0.39832 
##   b =  1.1112 +/- 0.10209 
## rho =  -0.39688

a) Vary the sampling resolution used to generate the data, keeping the same sampling range

x0      <- 0  # Signal peak
A.true  <- 2  # Signal amplitude
B.true  <- 1  # Background amplitude
Delta.t <- 5  # Exposure time

# signal widths
widths  <- c(0.1,  0.25,  1,  2, 3)

# sampling grid for computing posterior
alim <- c(0.0, 4.0)
blim <- c(0.5, 1.5)
Nsamp <- 100
uniGrid <- seq(from=1/(2*Nsamp), to=1-1/(2*Nsamp), by=1/Nsamp)

delta_a <- diff(alim)/Nsamp
delta_b <- diff(blim)/Nsamp
a <- alim[1] + diff(alim)*uniGrid
b <- blim[1] + diff(blim)*uniGrid

# initialize a list of signals
signals_w <- vector(mode = "list", length = length(widths))
sig.plot <- vector(mode = "list", length = length(widths))

for (i in seq_along(widths)){
  signals_w[[i]] <- posterior$new(seq(from=-7*widths[i], to=7*widths[i], by=0.5 *widths[i]), A.true, B.true, x0, widths[i], Delta.t)
  sig.plot[[i]] <- posterior$new(seq(from=-7*widths[i], to=7*widths[i], by=0.05*widths[i]), A.true, B.true, x0, widths[i], Delta.t)
}

# plot the signals, with the posteriors
i <- 1
while (i<=length(widths)){
  signals_w[[i]]$a.b.vectors(a, b)
  
  gg_title <- paste('w=', widths[i], sep='')
  gg_tmp <- (gg_signal(sig.plot[[i]], signals_w[[i]], gg_title)+signals_w[[i]]$plot.2d.posterior(gg_title))
  
  if (i==1) all_signals <- gg_tmp
  else all_signals <- all_signals / gg_tmp
  
  i <- i+1
}

all_signals+
  plot_layout(guides = "collect")+
  plot_annotation(
    title='Signal, background counts and unnormalized 2D posterior', 
    theme = theme(plot.title = element_text(size = 20, hjust = 0.5))
    )

# Compute marginalized posteriors, plot them, get parameters values
df_param <- data.frame('w'=widths)
params <- matrix(nrow = length(widths), ncol = 5)

i <- 1
while (i<=length(widths)){
  signals_w[[i]]$marginalization(delta_a, delta_b)
  
  gg_tmp <- (signals_w[[i]]$plot.marginalized(title=paste('w=', widths[i], sep='')))
  
  if (i==1) marg.distributions <- gg_tmp
  else marg.distributions <- marg.distributions / gg_tmp
  
  params[i,] <- signals_w[[i]]$param_values(delta_a, delta_b, print_param=FALSE, return_param=TRUE)

  i <- i+1
}

df_param['mean A'] <- params[, 1]
df_param['mean B'] <- params[, 2]
df_param['std A']  <- params[, 3]
df_param['std B']  <- params[, 4]
df_param['rho A, B'] <- params[, 5]

df_param %>% 
  kable(digits=4) %>%
  kable_classic(full_width=FALSE)

w	mean A	mean B	std A	std B	rho A, B
0.10	1.6302	1.1112	0.3983	0.1021	-0.3969
0.25	1.6302	1.1112	0.3983	0.1021	-0.3969
1.00	1.6302	1.1112	0.3983	0.1021	-0.3969
2.00	1.6302	1.1112	0.3983	0.1021	-0.3969
3.00	1.6302	1.1112	0.3983	0.1021	-0.3969

marg.distributions+
  plot_layout(guides = "collect")+
  plot_annotation(
    title='Marginalized posterior distributions', 
    theme = theme(plot.title = element_text(size = 20, hjust = 0.5))
    )

b) Change the ratio A/B used to simulate the data (keeping both positive in accordance with the prior)

In this point I keep \(B\) constant and equal to 1.

x0      <- 0  # Signal peak
w       <- 1  # Signal width
B.true  <- 1  # Background amplitude, I keep it constant
Delta.t <- 5  # Exposure time

xdat <- seq(from=-7*w, to=7*w, by=0.5 *w)
xplot <- seq(from=-7*w, to=7*w, by=0.05 *w)
  
# ratios A/B = A, because B=1
ratios  <- c(0.5, 1, 2, 4, 6)

# sampling grid for computing posterior
alim <- c(0.0, 6.0)
blim <- c(0.5, 1.5)
Nsamp <- 100
uniGrid <- seq(from=1/(2*Nsamp), to=1-1/(2*Nsamp), by=1/Nsamp)

delta_a <- diff(alim)/Nsamp
delta_b <- diff(blim)/Nsamp
a <- alim[1] + diff(alim)*uniGrid
b <- blim[1] + diff(blim)*uniGrid

# initialize a list of signals
signals_AB <- vector(mode = "list", length = length(ratios))
sig.plot <- vector(mode = "list", length = length(ratios))

for (i in seq_along(ratios)){
  signals_AB[[i]] <- posterior$new(xdat, ratios[i], B.true, x0, w, Delta.t)
  sig.plot[[i]] <- posterior$new(xplot, ratios[i], B.true, x0, w, Delta.t)
}

# plot the signals, with the posteriors
i <- 1
while (i<=length(ratios)){
  signals_AB[[i]]$a.b.vectors(a, b)
  
  gg_title <- paste('A/B=', ratios[i], sep='')
  gg_tmp <- (gg_signal(sig.plot[[i]], signals_AB[[i]], gg_title)+signals_AB[[i]]$plot.2d.posterior(gg_title))
  
  if (i==1) all_signals <- gg_tmp
  else all_signals <- all_signals / gg_tmp
  
  i <- i+1
}

all_signals+
  plot_layout(guides = "collect")+
  plot_annotation(
    title='Signal, background counts and unnormalized 2D posterior', 
    theme = theme(plot.title = element_text(size = 20, hjust = 0.5))
    )

# Compute marginalized posteriors, plot them, get parameters values
df_param <- data.frame('A/B'=ratios)
params <- matrix(nrow = length(ratios), ncol = 5)

i <- 1
while (i<=length(ratios)){
  signals_AB[[i]]$marginalization(delta_a, delta_b)
  
  gg_tmp <- (signals_AB[[i]]$plot.marginalized(title=paste('A/B=', ratios[i], sep='')))
  
  if (i==1) marg.distributions <- gg_tmp
  else marg.distributions <- marg.distributions / gg_tmp
  
  params[i,] <- signals_AB[[i]]$param_values(delta_a, delta_b, print_param=FALSE, return_param=TRUE)

  i <- i+1
}

df_param['mean A'] <- params[, 1]
df_param['mean B'] <- params[, 2]
df_param['std A']  <- params[, 3]
df_param['std B']  <- params[, 4]
df_param['rho A, B'] <- params[, 5]

df_param %>% 
  kable(digits=4) %>%
  kable_classic(full_width=FALSE)

A.B	mean A	mean B	std A	std B	rho A, B
0.5	0.2167	1.0798	0.1739	0.0900	-0.2391
1.0	0.5364	1.1072	0.2921	0.0976	-0.3872
2.0	1.6302	1.1112	0.3983	0.1021	-0.3969
4.0	3.5259	1.1557	0.4915	0.1041	-0.3375
6.0	5.2911	1.1662	0.4336	0.1038	-0.2563

marg.distributions+
  plot_layout(guides = "collect")+
  plot_annotation(
    title='Marginalized posterior distributions', 
    theme = theme(plot.title = element_text(size = 20, hjust = 0.5))
    )